#!/usr/bin/env python 

import sys
import os 
import logging

from Bio import SeqIO
from Bio import AlignIO
from Bio.SeqUtils import GC
from Bio.Align.Applications import ClustalwCommandline

def compute_GC(seqs_name, seqs): 
  gc_values = []
  for i in range(len(seqs)):
    gc_values.append(GC(seqs[i]))
  return gc_values

def transcription(seqs_name, seqs): 
  rna_seqs = [] 
  for i in range(len(seqs)):
    rna_seqs.append(seqs[i].transcribe())
  return rna_seqs

def translation(seqs_name, seqs): 
  prt_seqs = []
  for i in range(len(seqs)):
    prt_seqs.append(seqs[i].translate())
  return prt_seqs

def alligement(in_file_name, out_file_name, out_file_name_phy): 
  cline = ClustalwCommandline("clustalw2", infile = in_file_name, matrix = "BLOSUM", outfile = out_file_name)

  logging.info("Run clustalw")
  cline()
  logging.info("End clustalw")

  inp = open(out_file_name, "rU")
  out = open(out_file_name_phy, "w")

  alignments = AlignIO.parse(inp, "clustal")
  AlignIO.write(alignments, out, "phylip")

  inp.close()
  out.close()

def read_sequences(files): 
  seqs_name = []
  seqs = []
  for file_name in files: 
    type_file = ""
    if file_name.endswith("fasta"):
      type_file = "fasta"
    elif file_name.endswith("gb") or file_name.endswith("gbk"): 
      type_file = "genbank"

    count = 0

    for seq in SeqIO.parse(file_name, type_file):
      seqs_name.append(((file_name.split("/"))[::-1])[0])
      #seqs_name.append(seq.name)
      seqs.append(seq.seq)
      count += 1

    if count > 1: 
      logging.error("File %s contains more than one sequences", file_name)
      exit()
      
  return (seqs_name, seqs)  

if __name__ == "__main__": 
  logging.basicConfig(level=logging.INFO)    

  if len(sys.argv) <= 1:  
    sys.stderr.write("USAGE: main.py file1 file2 file3 ...\n")
    exit()

  prt_file = "temp.fasta"
  alg_file = "aligement.aln"
  phy_file = "aligement.phy"

  logging.info("Start doing job")
 
  logging.info("Start to checks format and exists files")

  for file_name in sys.argv[1:]:
    if not os.path.isfile(file_name):
      logging.error("File %s didn't exists", file_name) 
      exit()
    if not (file_name.endswith("fasta") or file_name.endswith("gb") or file_name.endswith("gbk")):
      logging.error("Bad format file %s.\nSupported follow format: file1.fasta or file1.gb or file1.gbk", file_name)
      exit()

  logging.info("Start to read sequences")
  seqs_name, seqs = read_sequences(sys.argv[1:])

  logging.info("Start to compute GC contents")
  gc_values = compute_GC(seqs_name, seqs)

  for i in range(len(seqs_name)):
    print "GC content " + seqs_name[i] + " = " + str(gc_values[i]) + "%"

  logging.info("Start to translate DNA to RNA")
  rna_seqs = transcription(seqs_name, seqs) 

  for i in range(len(seqs_name)):
    print "RNA for " + seqs_name[i] + ":\n" + rna_seqs[i]

  logging.info("Start to translate RNA to proteins")
  prt_seqs = translation(seqs_name, rna_seqs) 
  for i in range(len(seqs_name)):
    print "Protein in " + seqs_name[i] + ":\n" + prt_seqs[i]

  logging.info("Save proteins in temporary file %s", prt_file)
  for i in range(len(seqs_name)):
    with open(prt_file, 'a') as f:   
      f.write(">" + seqs_name[i] + "\n")
      f.write(str(prt_seqs[i]) + "\n")
   
  logging.info("Start to create aligement of proteins")
  alligement(prt_file, alg_file, phy_file)

  logging.info("Cleaning")
  os.remove(prt_file)
  logging.info("Finish")



